########文章用图
library(ggplot2)
library(dplyr)															######chromHMM染色质状态的p1
setwd("E:/5hmc_file/脑组织和GM12878的chromHMM数据")
anno=read.table("group意义.txt",head=T,sep="\t")
anno=data.frame(HMM.group=anno$STATE.NO.,description=anno$DESCRIPTION,num=anno$NUM,abbr=anno$MNEMONIC)

file=read.table("result_HMM_enrichment.csv",head=T,sep=",")
groupf=unique(file$con.file.name)

file1=file[file$con.file.name==groupf[2],]
file1$HMM.group=paste0("NO.",file1$HMM.group)
filee=merge(file1,anno,by="HMM.group")

col_names=c("HMM.group","OR","lower","upper","p.value")
result=data.frame(matrix(NA,15,ncol=5))
names(result)=col_names
group.HMM=unique(filee$HMM.group)
library(meta)
for(i in 1:15){
tdata=data.frame(HMM.group=filee$HMM.group,filee[,4:7])
tdata1=tdata[tdata$HMM.group==group.HMM[i],]
 tdata1$case_not=tdata1$case_in_region+tdata1$case_not
 tdata1$con_not=tdata1$con_in_region+tdata1$con_not
metaor3<-metabin(case_in_region,case_not,con_in_region,con_not,data=tdata1,sm="OR",studlab = HMM.group) 
OR=exp(metaor3$TE.fixed)
upper=exp(metaor3$upper.fixed)
lower=exp(metaor3$lower.fixed)
result[i,]$HMM.group =group.HMM[i]
result[i,]$OR=exp(metaor3$TE.fixed)
result[i,]$upper=exp(metaor3$upper.fixed)
result[i,]$lower=exp(metaor3$lower.fixed)
result[i,]$p.value=metaor3$pval.fixed
}
result=merge(result,anno,by="HMM.group")
 result$FDR=p.adjust(result$p.value,method = "bonferroni")
 result$group="nonsignificant"
 result[result$FDR<0.05,]$group="significant"
 result$OR=log2(result$OR)
 result$lower=log2(result$lower)
 result$upper=log2(result$upper)
 result=result[order(result$num),]
 result1=result[1:9,]
 result2=result[10:15,]
 result=rbind(arrange(result1,group,OR),arrange(result2,group,OR))
 result$num=1:15
 result1=result

p1=ggplot(data=result,aes(y=OR,x=num,group=group,color=group))+geom_errorbar(aes(ymin=lower,ymax=upper),width=0.1,color="black")+
    geom_point(aes(shape=group),size=3)+coord_flip()+geom_hline(yintercept = log2(1.5),color="red",linetype="dashed")+scale_x_continuous(breaks = result$num,labels = result$abbr)+
    theme_light()+ylab("log2(OR)")+geom_hline(yintercept = 0,color="black")+
    scale_color_manual(values = alpha(c("#3C5488","#E64B35"),0.9))+theme(panel.grid.minor = element_blank(),legend.position = "none")
	
													################组织特异性的图p2

setwd("E:/5hmc_file/组织特异性表达")
library(openxlsx)
file=read.xlsx("tissue_spec_enrichment_result_sig.xlsx",sheet=1)


result_c2=file[file$con.file=="nobias_AShM",]
 #result_c2=result_c2[result_c2$OR>1,]							######这部分是用于meta分析的
 #result_c2=result_c2[order(result_c2$OR),]
 #result_c2$case_not=result_c2$case_overlap+result_c2$case_not
 #result_c2$con_not=result_c2$con_overlap+result_c2$con_not
 dup_id=unique(result_c2[duplicated(result_c2$tissue),]$tissue)###把重复行挑出来重新处理求均值,然后求OR和pvalue
 result_c2_no_dup=result_c2[!result_c2$tissue %in% dup_id,]

 for(m in 1:length(dup_id)){
     rt_tmp =result_c2[result_c2$tissue==dup_id[m],]
     rt_tmp2=rt_tmp[1,]
     rt_tmp2$spec_filename=paste(unlist(gsub(".txt","",unique(rt_tmp$spec_filename))),collapse = "_")
     rt_tmp2$case_overlap=round(mean(rt_tmp$case_overlap))
     rt_tmp2$con_overlap=round(mean(rt_tmp$con_overlap))
     rt_tmp2[,8]=fisher.test(matrix(c(rt_tmp2[1,4],rt_tmp2[1,5],rt_tmp2[1,6],rt_tmp2[1,7]),nrow = 2))$estimate
     rt_tmp2[,9]=fisher.test(matrix(c(rt_tmp2[1,4],rt_tmp2[1,5],rt_tmp2[1,6],rt_tmp2[1,7]),nrow = 2))$p.value
     result_c2_no_dup=rbind(result_c2_no_dup,rt_tmp2)
 }
 result_c2_no_dup=result_c2_no_dup[order(result_c2_no_dup$OR),]
 result_c2_no_dup$FDR=p.adjust(result_c2_no_dup$p.value,method = "BH")
 result_c2_no_dup$tissue=gsub("[^[:alnum:]///' ]"," ",result_c2_no_dup$tissue)
 write.csv(result_c2_no_dup,"tis_sepc_enrichment_rt.csv",quote = F,row.names = F)
 rt=read.xlsx("tis_sepc_enrichment_rt.xlsx",sheet=1)
 rt$logpvalue=-log10(rt$p.value)
 rt$group="nosig"
 rt[rt$p.value<0.05,]$group="sig"
 p2=ggplot(rt,aes(y=OR,x=logpvalue,color=group))+geom_point(size=3)+scale_color_manual(values = alpha(c('#3C5488','#E64B35'),0.8))+
    labs(y="OR",x="-log10 (P-value)",title="tis sepc enrichment")+
    theme_classic(base_size = 12)+geom_vline(xintercept = -log10(0.05),color="red",linetype="dashed")+geom_hline(yintercept = 1,color="red",linetype="dashed")+
    geom_text_repel(data=subset(rt,rt$p.value<0.05),aes(label=tissname),size=3, fontface="bold",force = T,box.padding = unit(0.5, "lines"),point.padding = unit(0.8, "lines"), segment.color = "black", show.legend = FALSE)

 
 
 metaor_no_dup<-metabin(case_overlap,case_not,con_overlap,con_not,data=result_c2_no_dup,sm="OR",studlab = tissue) 
 forest(metaor_no_dup)
 p2_1=forest(metaor_no_dup)###森林图
 
 meta_data=as.data.frame(metaor_no_dup)
meta_data$issue=gsub("Brain...","",result_c2_no_dup$tissue)

meta_data=meta_data[order(meta_data$TE),]
meta_data$num=1:dim(meta_data)[1]
ggplot(data=meta_data,aes(y=TE,x=num,color=pval))+geom_point(size=3)+
    geom_errorbar(aes(ymin=lower,ymax=upper),width=0.1,color="black")+coord_flip()+scale_color_gradient(low = "red", high = "green")+
    scale_x_continuous(breaks = meta_data$num,labels = meta_data$issue)+theme_light()+ylab("log(OR)")+
    theme(panel.grid.minor = element_blank())				######erro bar还是太长了不好看，pass
	
															######有标记的散点图
library(ggrepel)
meta_data1=data.frame(issue=meta_data$issue,TE=meta_data$TE,num=meta_data$num)
 meta_data2=meta_data1
 meta_data2$TE=0
 meta_data_line=rbind(meta_data1,meta_data2)
 +geom_line(data = meta_data_line,aes(x=TE,y=num,group=issue))
p2_2=ggplot(meta_data,aes(TE,num))+geom_point(aes(size=3,color=pval))+geom_text_repel(aes(label = event.e))+
  scale_color_gradient(low = "#E64B35", high = "#3C5488")+scale_x_continuous(breaks=c(seq(0,2.5,0.5)),limits = c(0,2.5))+
  labs(x="log(OR)",y="issue")+scale_y_continuous(breaks = meta_data$num,labels = meta_data$issue)+
  theme_bw()+theme(panel.grid.minor = element_blank())
  
													################TWAS 富集，只挑精神类的p3
setwd("E:/0 公共数据库差异情况/db_for_5hmc/TWASdb/")
file=read.table("TWAS_enrichment.csv",head=T,sep=",")
file$FDR=p.adjust(file$p.value,method = "BH")
psy=c("AlzheimersProxyMetaIGAP_Marioni2018.dat","BD_Ruderfer2018.dat","BDSCZ_Ruderfer2018.dat","Depression_Nagel2018.dat","MDD_Wray2018.dat")
result3=file[file$TWAS.db %in%psy &file$con.file.name=="nobias_AShM.txt",]
result3$TWAS.db=gsub(".dat","",result3$TWAS.db)

anno=read.table("TWAS_status.txt",head=T,sep="\t")
result3=merge(anno,result3,by="TWAS.db")

metaor3<-metabin(case_in_region,case_not,con_in_region,con_not,data=result3,sm="OR",studlab = status) 
forest(metaor3)

meta_data=as.data.frame(metaor3)
meta_data=meta_data[order(meta_data$pval),]
meta_data$num=5:1
#meta_data$num=as.numeric(c(2,3,1,4,5))		#调整一下顺序
meta_data1=data.frame(status=meta_data$studlab,TE=meta_data$TE,num=meta_data$num)
meta_data2=meta_data1
meta_data2$TE=0
meta_data_line=rbind(meta_data1,meta_data2)

p3=ggplot(meta_data,aes(TE,num))+geom_point(aes(size=event.e,color=pval))+geom_text_repel(aes(label = event.e))+
    scale_color_gradient(low = "#E65B35", high = "#3C5488")+scale_x_continuous(breaks=c(seq(0,2,0.5)),limits = c(0,2))+
    labs(x="log(OR)",y="status")+scale_y_continuous(breaks = meta_data$num,limits=c(0,6),labels = meta_data$studlab)+
    theme_bw()+theme(panel.grid.minor = element_blank(),panel.grid.major = element_blank())
	 
	 
												###################TF富集图p4
setwd("E:/5hmc_file/motifbreakR_predict_enrich")
file=read.csv("enrichment_by_TF_predit_by_MBR.csv",head=T)
file2=file[file$FDR<0.1&file$OR>1,]
file2=file2[file2$con.file=="nobias_AShM",]
file2$`1/OR`=1/file2$OR
file2=file2[order(file2$OR,decreasing = T),]


file2$num=as.numeric(c(4:8,1:2,9:15,3,16:20))
file2=file2[order(file2$num,decreasing = T),]
file2=file2[order(file2$num),]
file2$num=20:1

data1=data.frame(TFname=file2$TFname,num=file2$num)####为画水平连线做准备
data1$`1/OR`=file2$`1/OR`
data2=data1
data2$`1/OR`=0
data_line=rbind(data1,data2)
+geom_line(data = data_line,aes(x=`1/OR`,y=num,group=TFname))
p4=ggplot(file2,aes(`1/OR`,num))+geom_point(aes(size=motif_AShM,color=FDR))+geom_text_repel(aes(label = motif_AShM))+
  scale_color_gradient(low = "#E64B35", high = "#3C5488")+labs(x="1/OR",y="TF.name")+
  scale_y_continuous(breaks = file2$num,labels = file2$TFname)+
  theme_bw()+theme(panel.grid.minor = element_blank())
  
  tmp=file2[1,]
  p5=ggplot(tmp,aes(`1/OR`,num))+geom_point(aes(size=motif_AShM))+scale_y_continuous(breaks=c(seq(19.5,20,0.5)),limits=c(19.5,20))+
  scale_x_continuous(breaks=c(seq(0,1,0.5)),limits = c(0,1))
  
  #####################################################TF富集图展示效果2
  file=read.csv("enrichment_by_TF_predit_by_MBR.csv",head=T)
file=file[file$con.file=="nobias_AShM",]
file[file$OR=="Inf",]$OR=27
library(ggplot2)
library(ggrepel)
file$group="nosig"
file[file$FDR<0.1,]$group="sig"
file$logFDR=-log10(file$FDR)
p4=ggplot(file,aes(y=OR,x=logFDR,color=group))+geom_point(size=3)+scale_color_manual(values = alpha(c('#3C5488','#E64B35'),0.7))+
    labs(y="OR",x="-log10 (FDR)",title="TF enrichment")+
    theme_classic(base_size = 12)+geom_vline(xintercept = -log10(0.1),color="red",linetype="dashed")+geom_hline(yintercept = 1,color="red",linetype="dashed")+
    geom_text_repel(data=subset(file,file$FDR<0.1),aes(label=TFname),size=3, fontface="bold",force = T,box.padding = unit(0.5, "lines"),point.padding = unit(0.8, "lines"), segment.color = "black", show.legend = FALSE)


  
#############################################################调控特征出现的可能性和羟甲基化升降一致性的统计
setwd("E:/5hmc_file/H3k的分析")
sel=list.files(pattern = "Brain")
library(ggplot2)
library(ggrepel)
library(RColorBrewer)

alias=c("AG","AC","GM","HM","ITL","MFL","SN","FBF","FBM")
for(i in 2:length(sel)){
ftp=read.csv(sel[i],head=T)
names(ftp)=c("term",names(ftp)[-c(1)])
ftp$group=alias[i]
file=rbind(file,ftp)
}
file$logpvalue=-log10(file$pvalue)

getPalette = colorRampPalette(brewer.pal(9, "Set1"))

p5=ggplot(file,aes(y=same_ratio,x=logpvalue,color=term))+geom_point(size=3)+scale_color_manual(values = alpha(c('#C71585','#D2691E','#E64B35','#3CB371','#4169E1','#708090','#000000'),0.9))+
  labs(y="same_ratio",x="-log10 (Pvalue)",title="H3k")+
  theme_classic(base_size = 12)+geom_vline(xintercept = -log10(0.05),color="red",linetype="dashed")+
  geom_text_repel(data=subset(file,file$pvalue<0.05),aes(label=group),size=3, fontface="bold",force = T,box.padding = unit(0.5, "lines"),point.padding = unit(0.8, "lines"), segment.color = "black", show.legend = FALSE)

  
  
  layout=matrix(c(1,1,2,2,1,1,2,2,1,1,3,3,1,1,3,3),4,4,byrow=T)
  
  layout=matrix(c(1,1,2,2,1,1,2,2,1,1,2,2,1,1,3,3,1,1,3,3,1,1,4,4,1,1,4,4,1,1,4,4,1,1,4,4),9,4,byrow=T)
  multiplot(plotlist=list(p1,p2,p4),layout=layout)
  
  multiplot(plotlist=list(p1,p5,p5,p5),layout=layout)
  
  multiplot(plotlist=list(p5,p2_2,p5,p4),layout=layout)
  
  multiplot(plotlist=list(p5,p5,p3,p5),layout=layout)